Conformal and covariant formulation of the Z4 system with constraint-violation damping 
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We present a new formulation of the Einstein equations based on a conformal and traceless decomposition of 
the covariant form of the Z4 system. This formulation combines the advantages of a conformal decomposition, 
such as the one used in the BSSNOK formulation (i.e. well-tested hyperbolic gauges, no need for excision, 
robustness to imperfect boundary conditions) with the advantages of a constraint-damped formulation, such as 
the generalized harmonic one (i.e. exponential decay of constraint violations when these are produced). We 
validate the new set of equations through standard tests and by evolving binary black hole systems. Overall, the 
new conformal formulation leads to a better behavior of the constraint equations and a rapid suppression of the 
violations when they occur. The changes necessary to implement the new conformal formulation in standard 
BSSNOK codes are very small as are the additional computational costs. 

PACS numbers: 04.25.D-, 04.25.dg 



I. INTRODUCTION 

Numerical relativity has seen, over the last few years, a 
truly remarkable development. Starting from the first simu- 
lations showingthat black-hole binaries could be evolved for 
a few orbits [QJ-yl], or that black holes could be produced from 
unstable stellar configurations using simple gauges and with- 
out excision fl], new results have been obtained steadily. As a 
result, it is now possible to simulate binary black holes |5[] and 
binary neutron stars 16|] accurately for dozens of orbits, from 
the weak-field inspiral, down to the final black-hole ringdown 
(see also ^ for recent reviews on binary black holes and 
neutron stars, respectively). In addition, the progress in nu- 
merical relativity has also been accompanied by a comparable 
progress of analytical approximation techniques, which have 
been shown to be able to reproduce the numerical results to 
very high precision both for binary black holes J^, [T^] and 
for binary neutron stars ifllll . Finally, numerical simulations 
have now investigated scenarios never considered before and 
that could lead to a new and deeper understanding of the as- 
trophysics of compact objects |L?[ [LiTl . 

There are several reasons behind this rapid progress, and 
the use of more accurate numerical techniques and the avail- 
ability of larger computational facilities are certainly among 
the most important ones. None of these, however, would be 
useful without the use of formulations of the Einstein equa- 
tions that are well-suited for numerical evolutions. Most of the 
present three-dimensional (3D) numerical-relativity codes im- 
plement either one of the two formulations discussed below. 
The first and most popular one is the conformal and traceless 
reformulation of the 3 + 1 ADM equations Il4tl . which is also 
known as the BSSNOK (or BSSN) formulation [SHI- The 
second formulation is instead based on the use of a fully 4D 
form of the Einstein equations in coordinates that resemble 
the harmonic ones and is therefore known as the Generalized- 
Harmonic formulation (GH) lfl8ll . 



There are several differences between these two formu- 
lations, each having its own advantages and disadvantages. 
One of the main advantages of BSSNOK is that, being based 
on a conformal decomposition, it can separate potential sin- 
gular terms in the conformal factor. In addition, it can 
count on well-tested and robust gauge conditions, such as the 
singularity-avoiding slicing conditions of the 1 + log fam- 
ily lHH]. Similarly, the spatial gauges can rely on the hy- 
perbolic Gamma-driver condition for the shift vector luOll 
(or some recent variants for unequal-mass binaries I2H - I23ID . 
which removes to a large extent, the gauge dynamics near the 
compact objects. When combined, these two gauge choices 
eliminate the need to excise a region of the computation do- 
main inside the apparent horizon, greatly simplifying the nu- 
merical infrastructure. Finally, the use of the momentum con- 
straint equations (but not of the energy constraint) in the evo- 
lution of the dynamical variables, which is crucial for ensur- 
ing strong hyperbolicity, provides BSSNOK with a certain 
"forgiveness", so that the violation of the constraints does 
not grow rapidly, even when boundary conditions which are 
constraint-violating are used near the strong-field region. 

In contrast, the GH formulation uses a generalized har- 
monic gauge which cannot deal with the physical singularity 
inside the apparent horizon. As a result, at least for the gauges 
considered so far (see also ll24ll25ll ), it requires the use of ex- 
cision and thus of numerical techniques that are devised for 
handling a special region of the computational domain l26ll . 
To its advantage, however, the GH formulation leads to a set of 
equations whose principal parts are wave equations and thus 
with very well-known mathematical properties. In addition, 
the use of damping terms allows for the dynamical control of 
the constraint violations and thus for a powerful way of reduc- 
ing them when necessary. Of course, a solution with smaller 
constraint violations will intrinsically be a more accurate so- 
lution to the Einstein equations. 

Clearly, it would be useful to employ a formulation of the 
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Einstein equations that combines the best of both worlds and 
thus that has the robustness and gauge conditions of the BSS- 
NOK formulation but, at the same time, has well-defined 
mathematical properties and the possibility of dynamically 
controlling the constraint violations as the GH formulation. 
As we will show, these properties are met by a new conformal 
and covariant formulation of the Z4 system with constraint- 
violation damping. This is obtained by starting from the fully 
covariant Z4 formulation (27| and by performing a confor- 
mal decomposition which includes all the nonprincipal terms 
coming from the covariant form of the equations. In addition, 
damping terms are included for controlling the constraints 
in the spirit of the GH formulation. We will refer to this 
new formulation as the conformal and covariant Z4 system, 
i.e. CCZ4, and present tests of its behavior by considering 
evolutions in vacuum of gauge waves in ID and isolated and 
binary black holes in 3D. 

It should be remarked that this is not the first time that 
a conformal decomposition of the Z4 system has been pro- 
posed and indeed a very interesting attempt has been made in 
Ref. lf28tl . where it was named Z4c. Although the tests pre- 
sented in Ref. l28ll were performed in spherical symmetry, 
they already highlighted the potential of a conformal formu- 
lation of the Z4 system, especially in the presence of matter 
(see also l29l l30lo . Unfortunately, we were not able to ob- 
tain equally good results when evolving the formulation of 
Ref. 12811 in vacuum and in 3D; at the same time, we did not 
find that our CCZ4 formulation is more sensitive to boundary 
problems than the BSSNOK one (this was a point raised in 
Ref. JH). 

The structure of the paper is as follows. In Sec. [Til we de- 
rive the full set of the CCZ4 equations starting from the co- 
variant form of the Z4 system. In Sec. [HI] we introduce the 
details of the numerical infrastructure and present a numeri- 
cal comparison between the CCZ4 and the BSSNOK systems 
for a gauge-wave test and for binary black-hole simulations. 
Finally, the conclusions are summarized in Sec.lIVI 

II. THE CONFORMAL COVARIANT Z4 SYSTEM 

The Z4 formulation was introduced as a covariant exten- 
sion of the Einstein equations 12711 . where the original elliptic 
constraints are converted into algebraic conditions for a new 
four- vector Z M . This formulation can be derived from the co- 
variant Lagrangian 

C = g»"[R lu , + 2V M Z v ], (1) 



by means of a Palatini-type variational principle l3lil . The 
vector Z^ measures the deviation from the Einstein field equa- 
tions. The algebraic constraints Z^ = amount therefore to 
the fulfilling of the standard energy-momentum constraints. 
In order to control these constraints, the original system was 
supplemented with damping terms such that the true Einstein 
solutions (i.e. the ones satisfying the constraints) become an 
attractor of the enlarged set of solutions of the Z4 system l32ll . 
The Z4 damped formalism can be written in covariant form as 

R 

- (1 + K 2 )g^n a Z"} = 8tt(T^ - fa^T) , (2) 

where n M is the unit normal to the time slicing, the stress- 
energy tensor and T its trace, i.e. T = g^T^ . The (constant) 
coefficients m are free parameters related to the characteristic 
time of the exponential damping of constraint violations. As- 
suming energy-momentum tensor conservation, the Bianchi 
identities lead to the constraint-propagation system 

V^V ' yZ^ + R^Z" = -KiV y [n tl Z u + n„Z tl + K 2 g p .un a Z' J ] . 

(3) 

It has been shown in Ref. Il32ll that all the constraint-related 
modes are damped when 

Ki > k 2 > -1. (4) 

The Z4 formulation can be rewritten as a Cauchy problem 
by performing the 3 + 1 decomposition of the spacetime, in 
which the line element reads 

ds 2 = -a 2 dt 2 + 7 y (dx l + ftdt) (dx j + ftdt) , (5) 

where a is the lapse function, /3 1 is the shift vector and 7^ the 
intrinsic metric of the constant-time slices. The Einstein equa- 
tions within this decomposition lead to the well-known ADM 
system 03], which is usually cast as a system of evolution 
equations for the extrinsic curvature Kij and the three-metric 
7ij, plus four elliptic equations for the energy (or Hamilto- 
nian) and the momentum constraints, involving space deriva- 
tives of the dynamical fields 7^ and K%j. In the Z4 for- 
mulation, the energy-momentum constraints become evolu- 
tion equations for Z^, modifying the principal part of the 
ADM system and converting it from weakly to strongly hy- 
perbolic 1I33I1 . The 3 + 1 decomposition of the Z4 formulation 
including the damping terms reads 



{dt 



-2aKi 



-8ira 



R,j + ViZj + WjZi - 2 Ki Kij + (K - 26) K tj - Ki (1 + k 2 )Q Jij 



Sij --(3-t) 7y 



(6) 



(7) 



3 



(d t - l p ) e = 5 



A + 2 Vj^ + (A' - 29) A' - K ij Kij - 2 



2«l(2 + K 2 ) 6 - 167TT 



- = a [ Vj (A; 3 - <V A") + d& - 2 A/ Z 3 -Q-±- «iZ< - 8tt # ] , 



(8) 
(9) 



where Cp is the Lie derivative along the shift vector /3, is the 
projection of the Z4 four-vector along the normal direction, 
= n^Z^ = aZ°, and the following definitions apply for 



we express the metric 7,^ in terms of a conformal metric 



1 ij 



7ij with unit determinant <fr = (det(7jj)) 



-1/6 



while the extrinsic curvature Kij is decomposed into its trace 



A, 



4> 2 (K ij --K lij ). 



matter-related quantities r = n^nJT^ ', S t = n v T\, S tJ = A = Kijj ij and in its trace-free components 
T 

Equations ©-(O must be complemented with suitable 
gauge conditions that determine the system of coordinates 
used during the evolution. Of all the possible options, the 
most interesting ones are those which preserve the hyperbol- 
icity of the full evolution system, such as the 1 + log family 
and the Gamma-driver shift condition. 

As a first step towards deriving the CCZ4 formulation, 



(10) 



This allows us to write the three-dimensional Ricci tensor as 



Ft 



Ri 



ij t R-ij' tnus splitting it into a part containing con- 
formal terms and another one containing space derivatives of 
the conformal metric 



z,li 



Rij = -ifi^did^j + 7fe(^) r + r T (m + 7 



Rj 



(11) 
(12) 



r 



where 



(13) 



The conformal and covariant Z4 formulation (CCZ4) is thus 
given by the following system of evolution equations 



d t % = -2aA™ + 2%^ (3 k -pah P k + fi k d k % , (14) 
d t A l3 = 2 [-ViVj-a + a (A„- + V t Zj + VjZi - 87r5y)] TP + a^ij ( K ~ 2e ) 

-2aA a A] + 2A Hi d j} f3 k - ^A io d k (3 k + p k 8 k A l3 , (15) 

d t <t> = p^K - ^4>d k p k + (3 k d k , (16) 

d t K = -V i V l a + a(R + 2V t Z l + A 2 - 26A) + ftdjK - 3aKi (1 + k 2 ) 9 + Aira (S - 3t) , (17) 

d t O = (r+ 2\7 l Z t - A VJ A li + ^A 2 - 26A^j - Z%a + /3 k d k Q - a«i (2 + k 2 ) 9 - 87rar , (18) 

<9 t f * = 2a (f) k A^ - 3A ij ^ - p^djKj + 2 7 fci (ad k 6 - <dd k a - %aKZ k ^ - 2A ii d j a 

+i kl d k d l p + p ik d k d t p l + h l d k p k - f k d k p + 2k 3 Q-"z ; <x - i ]k z 3 d k p' 

+P k d k r - 2aK l f ' } Z, J - Vonaf* Sj , (19) 

dtot = -2a (K - 29) + f3 k d k a , (20) 

d t p = fB l + (3 k d k p\ (21) 

dtB 1 = d t f l -l3 k d k r +f3 k d k B l -r 1 B\ (22) 



where we have defined Note that the choice made with the definition d23l ) is equiva- 

r = r + 2f j Zj . (23) 
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lent, in the ADM context, to adding the momentum constraint 
to the right-hand-side of the evolution equation of P. In the 
context of the Z4 formulation, this just amounts to replacing 
the vector Zj by the quantities P in the set of basic fields to 
be evolved. 

The gauge conditions (120b — (f22l> correspond respectively to 
the standard "1 + log" slicing condition and to the original 
form of the gamma-driver shift condition, where a generic 
gauge parameter / was introduced |20]. Note that in the Z4 
formulation there is an additional propagation speed and the 
standard BSSNOK choice of / = 3/4 can then lead to weak 
hyperbolicity when the lapse a is close to 1 . This is why safer 
choices, such as / = 1, have been proposed in Ref. 112811 . In 
this paper we use / = 3/4 to be as close as possible to a 
standard BSSNOK formulation, but we also consider how the 
system of equations reacts when switching to / = 1. 

We also note that experimentation with black-hole space- 
times and the emergence of unstable behaviors, has induced us 
to introduce an extra parameter, « 3 , affecting some quadratic 
terms in the evolution Eq. ( fT9] l for T l . As discussed before, 
this equation corresponds to the evolution of Zi, so this is 
not just a gauge choice, but rather an essential ingredient of 
the Z4 system. Indeed, the covariance inherent to the confor- 
mal decomposition of the Z4 system is broken unless we take 
Ks = 1. For some of the tests presented in this paper we re- 
tain a fully covariant formulation (i.e. with K3 = 1). However, 
this is not possible for black-hole spacetimes, where nonlinear 
couplings with the damping terms, which are important for 
reducing the violations in the constraints, lead to numerical 
instabilities. As a result, for black-hole spacetimes we have 
resorted to a noncovariant and conformal formulation of the 
Z4 system (i.e. with K3 = 1/2) (see discussion in Sec. IIIIBl 
for details). 

A number of remarks are important at this point. First, al- 
though the structure of the CCZ4 formulation is very simi- 
lar to the BSSNOK one, there is an important difference in 
the evolution of the trace-free variable Aij. In the BSSNOK 
formulation, in fact, the Hamiltonian constraint is assumed 
to be satisfied exactly and thus used to eliminate the Ricci 
scalar from the right-hand-side of the evolution equation for 
Aij lf20ll . In the CCZ4 system, on the other hand, the evolution 
of follows directly from (the trace-free part of) the orig- 
inal ADM evolution equation for the extrinsic curvature Kij, 
plus the extra terms in Zi and 0. Second, the equivalent of 
the trace of the extrinsic curvature in BSSNOK formulations 
is given by 



K 



K - 2 6 , 



(24) 



again because the Hamiltonian constraint is assumed to re- 
move the Ricci scalar from the evolution equations in the 
BSSNOK approach. In the CCZ4 system, we rather use (the 
trace part of) the ADM evolution equation for Kij, modulo 
some Zi and terms. 

A closer look at the resulting CCZ4 system shows that it 
is not fully equivalent to the Z4 system, modulo a rearrange- 
ment of the dynamical fields. There are two extra fields which 
were not present in the Z4 system, namely dot 7^ and tr A^ . 



These are not dynamical fields at the continuum level, where 
the consistency constraints 



dot 7y = 1 



tr At 



0. 



(25) 



hold by construction. But at the discrete level, these are just 
two more constraints, which can be dealt with in many differ- 
ent ways. For instance: 

• Constrained approach. We could enforce d25l l at ev- 
ery integration step, by removing the trace of Aij 
and rescaling 7^- as it is usually done in BSSNOK 
codes [34]. The remaining dynamical modes have then 
the same characteristic structure of the original Z4 sys- 
tem. This is the safest choice, and we will use it in the 
tests presented in this paper. 

• Relaxed approach. We could instead relax ( |25] l, en- 
forcing it just on the initial/boundary data. In this way 
the two extra dynamical modes propagate along nor- 
mal lines, as their evolution equations [i.e. the trace of 
Eqs. (fT4l-([T5|>l are trivial. Note that in this case the 
trace of the first term in the evolution Eq. (fT~4-b must 
be removed explicitly to avoid any spurious numerical 
modes by evolving: 



<9t7i? 



-2a A. 



1 



-kl 



Moreover, in tests like the robust stability or the gauge 
waves, it may be necessary to keep also under control 
the trace of Aij . This can be achieved by adding, for 
instance, a damping term proportional to 7^ trAy to 
the evolution Eq. (15[ . 



Finally, the ADM constraints are given by 



H 



R - Ki.K ij + K 2 



(26) 



Mi = -f\diKij - diKji - TjiK mi + r™K m i) . (27) 

In the results presented below we compute the constraint vi- 
olations for both the BSSNOK and CCZ4 systems using the 
ADM quantities computed from the evolution variables corre- 
sponding to the two systems, allowing for the correspondence 
(El). 



III. NUMERICAL RESULTS 

In this section we validate the robustness and accuracy of 
the CCZ4 evolution system and compare it against the BSS- 
NOK system in two different cases: the gauge-waves test and 
black-hole spacetimes. In addition, we have performed sev- 
eral evolutions with the robust-stability test to ensure that the 
system is stable to linear perturbations, recovering the ex- 
pected results (see l35ll for a discussion of this test). 

The numerical setup used in the simulations presented 
here is the same one discussed in Ref. l36Tl and more re- 
cently applied to the Llama code described in Ref. l37ll . 
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The latter makes use of higher-order finite-difference algo- 
rithms satisfying the summation-by-parts rule (up to 8th or- 
der in space) and a multiblock structure for the outer com- 
putational domain. More specifically, we use a central cu- 
bical Cartesian patch containing multiple levels of adaptive 
mesh refinement with higher-resolution boxes. The Carte- 
sian grid is surrounded by 6 additional patches with the 
grid points arranged in a spherical-type geometry, with con- 
stant angular resolution to best match the resolution require- 
ments of radially outgoing waves. This allows us to move 
the outer boundary to a radius where it is causally discon- 
nected from the binary at a tiny fraction of the computational 
cost which would be necessary to achieve the same resolu- 
tion with a purely Cartesian code. The time evolution is 
based on the method-of-lines with a 4th order Runge-Kutta 
algorithm. Our general computational infrastructure is based 
on the Cactus framework and we are using pa ckag es such 
as TwoPunctures IB8[|. AHFinderDirect 113911 and of 
SummationByParts 14011 . which are freely available and 
part of the Einstein Toolkit. In addition, our evolutions make 
use of the mesh-refinement driver Carpet ll4lll , which imple- 
ments higher-resolution boxes with multiple levels of adaptive 
mesh refinement. 



A. Gauge Waves 

A classical test for different formulations of the Einstein 
equations is offered by the "gauge-wave" IT35tl . in which 
a fictitious one-dimensional pulse propagating along the in- 
direction can be simulated by performing a conformal trans- 
formation of the Minkowski metric in the two-dimensional 
sector spanned by the (t, x) coordinates, namely using the line 
element 



ds 1 



h(x, t) (~dt 2 + dx 2 ) + dy 2 + dz 1 



(28) 



The solution of the pulse at any time is just given by the 
advection of the initial profile of the gauge wave, which can 
be set to be smooth and periodic by choosing a sine-like initial 
data of the type jUl 



h(x,t = 0) = 1 - A sin 



(29) 



with an amplitude A < 1. Although this test is apparently triv- 
ial as it does not involve the solution of the Einstein equations 
in a very nonlinear regime, it nevertheless represents a seri- 
ous benchmark even for formulations as robust as BSSNOK, 
which indeed does not pass it ll42tl . 

Following i42tl . we choose an amplitude of A = 0.1 in 
a domain of L = 1 with three uniform resolutions ho/L = 
{1/50, 1/100, 1/200} and periodic boundary conditions. No- 
tice that the metric form d28l ) corresponds to an harmonic 
slicing condition with zero shift, so we have to change our 
preferred coordinate choice (i.e. the 1 + log slicing with the 
Gamma-driver) to perform this test. Furthermore, different 
implementations of the CCZ4 formulation: one in which 
the constraints are damped with coefficients m — 1/L and 
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FIG. 1 : L-infinity norm of the Hamiltonian constraint in the gauge- 
wave test, when performed with a CCZ4 formulation with damping 
terms (black solid line), with a CCZ4 formulation without damp- 
ing terms (blue dotted line), or with the BSSNOK formulation (red 
dashed line). Clearly, the Z4u and the BSSNOK formulations are 
unstable (cf. Fig. 5 of Ref. ll42ll ) and a similar behavior will be en- 
countered also in black-hole spacetimes (cf. Figure^. 



K2 = 0, and one in which the constraints are undamped, 
i.e. K\ = = K2- We will refer to these two cases as to "Z4d" 
and "Z4u", respectively (Note that in these tests the shift is set 
to zero and hence we do not need to specify a value for Ks, 
which we take to be one). 

The infinity-norm of the Hamiltonian constraint relative to 
simulations at the highest resolution is displayed in Fig. Q] 
for the damped CCZ4 formulation (black solid line), for the 
undamped CCZ4 formulation (blue dashed line), and for the 
BSSNOK formulation (red dotted line). Clearly, the BSSNOK 
and the CCZ4 formulation without damping terms fail before 
50 crossing times (BSSNOK after 42 crossing times and Z4u 
after 56 crossing times) as indicated by the an exponential 
increase in the violation of the energy constraint. However, 
with the addition of the damping terms, the CCZ4 formula- 
tion is able to accurately evolve this solution for more than 
1000 crossing times, while preserving the profile of the pulse. 
Furthermore, we have verified that the evolved solution con- 
verges to the expected spatial-discretization order (i.e. either 
4th or 8th order), with only a very small phase error when 
using the 8th order scheme. 

Overall, this test shows that the dynamical control of the 
energy constraint via the damping term K\ is crucial to attain 
a stable evolution, even in such a simple type of spacetimes. 
We also note that this test is more demanding for conformal 
formulations, where there is more than one component of the 
metric which is nontrivial. This is confirmed by comparing 
our results with those in Ref. j43ll . where the standard Z4 for- 
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TABLE I: Properties of the black-hole binaries simulated. The first column indicates the type of outer boundary and whether causally con- 
nected, ho is the grid spacing on the coarsest Cartesian grid, which is equal in all cases to the radial grid spacing in the angular patches. 
A^ng is the number of cells in the angular directions in the angular patches. Ri n and R ou t are the inner and outer radii of the angular patches. 
TVicv. is the number of refinement levels (including the coarsest) on the Cartesian grid, and 2 n cv indicates the size of the cubical refinement 
boxes centered on each black hole. The unit of the spacetime mass M is chosen such that each black hole has mass 0.5M in both the single 
and binary black cases. Finally, the last three columns contain the relative difference in the £ — m = 2 gravitational- wave phase between 
evolutions carried out with either the BSSNOK formulation (0b), the CCZ4 formulation with damping terms (0z4d), or the CCZ4 formulation 
without damping terms (0z4u). 



mulation, i.e. not implementing a conformal decomposition, 
was able to pass this test without the need of damping terms. 
The GH formulation also passes this test. 



B. Black-Hole Spacetimes 

Before considering black-hole binaries, we have tested ex- 
tensively our new CCZ4 formulation in the evolution of single 
nonspinning black holes. This has allowed us to determine 
how different choices for the damping coefficients m and K2 
influence the solution and, in particular, the violation of both 
the ADM and the constraints. In this way we have con- 
cluded that most of the dynamics in the evolution of the con- 
straint equations comes from the first damping coefficient, so 
that K2 — represents a sensible choice and is the one that we 
will consider hereafter. On the other hand, increasing values 
of Ki produce lower violations of the constraints and a value 
of Ki « 0.1/M seems optimal in this sense. Higher values, in 
fact, lead only to marginal improvements of the solution, but 
also tend to increase the stiffness of the damping terms. 

An important and unexpected result obtained when imple- 
menting the CCZ4 formulation in black-hole spacetimes is 
that subtle and nonlinear couplings can occur, leading to un- 
stable evolutions also for those choices of the coefficients that 
are perfectly stable in other spacetimes. While, in fact, we 
have carried out stable evolutions of the robust-stability test 
with the covariant and damped CCZ4 formulation (i.e. with 
K3 = 1 and Ki 7^ 0), we were not able to obtain stable evo- 
lutions of black-hole spacetimes with K3 = 1, although the 
growth time of the instability does change with the values of 
K\ (see discussion around Fig. |4). Clearly, nontrivial cou- 
plings seem to appear between these coefficients, which de- 
pend on the degree of nonlinearity and which deserve further 
investigation to be properly understood. 

On the whole, and as we will detail below, we have 
found that accurate and stable evolutions of binary black-hole 
spacetimes can be obtained with the damped noncovariant Z4 
systems (i.e. with Ks = 1/2, K\ = 0.1/M). On the other 
hand, covariant and conformal Z4 formulations that are either 



damped (i.e. with K3 = 1, K\ ^ 0), or undamped (i.e. with 
K3 = 1, K\ = 0), have been found to lead to unstable evolu- 
tions, although on rather different timescales and with variable 
degree of accuracy (see discussion below). 

The initial data of the binary black-hole evolutions is 
obtained from a circular-orbit condition at the third post- 
Newtonian order PPill and corresponds to an equal-mass non- 
spinning binary with an initial coordinate separation of D = 
8 M. The binary performs about 3.5 orbits before merging 
and settles to an isolated spinning black hole after t rj 360 M. 
To carry out a meaningful comparison, the binary is evolved 
with the BSSNOK and the CCZ4 formulations keeping the 
same choice for the gauges, namely the 1 + log slicing con- 
dition and the Gamma-driver shift condition with / = 3/4, 
77 = 2/M, and the same grid setup. For the latter, in par- 
ticular, we have considered three different choices aimed at 
determining the influence of the outer boundaries on the qual- 
ity of the solution. This is a point discussed in Refs. ll28[|29ll . 
where it was argued that the Z4c formulation is more sensitive 
than the BSSNOK one to incorrect (or constraint-violating) 
boundary conditions. As a result, we consider three different 
classes of simulations depending on the treatment of the outer 
boundary: (i) multiblock padding and spherical outer bound- 
ary which is causally disconnected (i.e. at ~ 2200 M for 
a simulation lasting ~ 800 M); (ii) multiblock padding and 
spherical outer boundary which is causally connected (i.e. at 
~ 350 M); (iii) Cartesian outer boundary which is causally 
connected (i.e. at ~ 200 M). For case (i), we reduce the or- 
der of the finite-difference operator at the outer boundary but, 
because it is causally disconnected, the initial conditions are 
preserved there. For case (ii), instead, we impose reflecting 
boundary conditions so as to "stress" the solution with data 
from the outer boundary which is constraint-violating and in- 
jected mostly at the time of the reflection. Finally, in case 
(iii) we have applied ordinary, outgoing Sommerfeld bound- 
ary conditions to all variables, again triggering violations in 
the constraint equations. 

All the properties of the grid structure and the treatment of 
the outer boundary are summarized in Table U] where ho is 
the grid spacing on the coarsest Cartesian grid, which is equal 
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FIG. 2: Real part of the £ = m = 2 mode of the gravitational wave- 
form *4 for an equal-mass nonspinning black-hole binary. Differ- 
ent lines refer to evolutions with the noncovariant formulation with 
and without damping terms, i.e. with K3 = 1/2 and Ki = 0.1/M, 
K2 = (Z4d), or K3 = 1/2 and n\ — «2 = (Z4u). The two evo- 
lutions are indicated, respectively, as Z4d and with a black solid line 
or as Z4u and with a blue dotted line; the BSSNOK formulation is 
shown with a red dashed line. Shown in the inset is a magnification 
of the merger. 



in all cases to the radial grid spacing in the angular patches. 
./Vang is the number of cells in the angular directions in the an- 
gular patches, while i?i n and R out are the inner and outer radii 
of the angular patches, respectively. In the case of a Carte- 
sian outer boundary, i? out represents the distance to the outer 
boundary along coordinate lines. Finally, iVie V , is the number 
of refinement levels (including the coarsest) on the Cartesian 
grid, while 2 n ev indicates the size of the cubical refinement 
boxes centered on each black hole. 

As final remark before discussing the results, we note that 
all the rest being the same, at any given resolution the CCZ4 
system has a smaller violation of the constraints than the BSS- 
NOK one. At the same time, however, because the violations 
of both the energy and momentum constraints are part of the 
evolution equations in the CCZ4 system, the latter is more 
strongly affected than BSSNOK one, for which only the vio- 
lations of the momentum constraint are included in the evo- 
lution system. As a result, the CCZ4 formulation requires a 
comparatively higher minimum-resolution treshold in order to 
enter a convergent regime. 

A first comparison of the behavior of the different formula- 
tions is offered in Fig. [2] where we show the I = m = 2 mode 
of the gravitational waveform $4 as extracted on a sphere of 
coordinate radius r = 100 M (see ll37ll for details on the ex- 
traction procedure). Different lines refer to simulations us- 
ing either the noncovariant formulation with damping terms, 



FIG. 3: Differences in the phase evolutions at the high, medium and 
low resolutions, respectively (these are indicated as HR, MR and 
LR). The top panel refers to the BSSNOK formulation, while the bot- 
tom one the the noncovariant damped CCZ4 formulation (Z4d). The 
differences between the low and medium resolutions are also scaled 
with the appropriate convergence coefficients (marked as CF4 and 
CF8, see text) to highlight the convergence order of the solution; all 
the data refers to simulations with a multiblock padding and causally 
disconnected outer boundary. Note that at these resolutions the CCZ4 
formulation has larger phase errors, but due its higher convergence 
factor, these errors are expected to decay at a faster rate than for 
BSSNOK. 



i.e. with k 3 = 1/2 and Ki = 0.1/M, k 2 = (Z4d, black 
solid line), or to the noncovariant formulation without damp- 
ing terms, i.e. with K3 = 1/2 and K\ = K2 — (Z4u, blue 
dotted line). Also shown as a reference is a simulation with 
the BSSNOK formulation (red dashed line) using the same 
numerical setup. The simulations refer to the highest resolu- 
tion (i.e. ho/M = 0.48) and the grid having the multiblock 
padding and an outer boundary at i? ou t = 2192.16 M, 

The first obvious thing to note is that all simulations lead 
to a stable merger and ringdown at all the resolutions consid- 
ered. Furthermore, while a small phase difference is present 
between the Z4 and the BSSNOK runs, this difference is very 
small and Acj> < 0.02 rad over the whole simulation. As a 
comparison, the phase difference between the Z4 and the Z4u 
simulations is Acj) < 0.002 rad (see Table [I] for the relative 
maximum differences). 

Although the phase differences between the waveforms ob- 
tained with the two formulations is relatively small, it also 
decreases with the resolution, thus indicating that both for- 
mulations would yield the same phase evolution in the con- 
tinuum limit. The rate of convergence, however, is different 
when considering either the BSSNOK or the CCZ4 formula- 
tion. This is shown in Fig. [3j where we report the residuals 
in the phase evolutions at the high, medium and low resolu- 
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tions, respectively (these are indicated as "HR", "MR" and 
"LR"). The differences between the low and medium resolu- 
tions are also scaled to highlight the convergence order of the 
solution. More specifically, the HR, MR and LR refer to simu- 
lations with the coarsest resolutions of ho/M = 0.6, 0.48, 0.4 
(cf. Table U). The convergence coefficients corresponding to 
these resolutions and used for rescaling are CF4 = 3.0898, 
for a convergence factor of 4.5 in the BSSNOK case, and 
CF8 = 7. 1906 for a convergence factor of 8.5 in the Z4d case. 
Note however that, as mentioned above, the CCZ4 formula- 
tion needs a higher resolution to enter the convergence regime, 
while a triplet of resolutions with ho/M = 0.8,0.6,0.48 
would be enough to show convergence at about 4th order for 
the BSSNOK runs. 

Beside this minimum resolution threshold, the additional 
computational expenses required by the CCZ4 formulations 
are not significant. The difference with the BSSNOK sys- 
tem consists in an additional evolution equation for the scalar 
variable 0, which would amount to solving 25 evolution equa- 
tions (instead of 24 as in BSSN), implying around 4% higher 
computational costs. However this is an over-estimate, as in 
reality the time spent in computing the evolution equations 
depends on the computational infrastructure. In our case, it 
is about half of the total time of a binary black-hole simu- 
lation, while the other half is dedicated to mesh-refinement, 
gravitational-wave extraction and other analysis routines. 

All in all, we find that for the highest resolutions used the 
results of the BSSNOK runs converge at about 4th order (top 
panel in Fig. [3}^ while the Z4d runs converge at about 8th or- 
der (bottom panel in Fig. [3]); in both cases, the convergence or- 
der is lost in the very final stages of the merger. It is a present 
unclear why the two formulations yield, with the same com- 
putational infrastructure, two different convergence rates. It is 
possible that the constraint-damping properties of the CCZ4 
formulation are able to suppress the small violations coming 
from the reflections across refinement boundaries, that are a 
major source of error and one of the largest obstacles to attain 
clean convergence. However, more efforts (and considerable 
computational costs) need to be invested to assess whether this 
is the correct explanation. 

A useful way to appreciate the different behavior of the 
two formulations is shown in Fig. [4] which reports the evo- 
lution of the L2-norm of the ADM energy (i.e. the violation 
of the Hamiltonian constraint) for the covariant CCZ4 formu- 
lation with and without damping (light-blue dot-dashed line 
and magenta long-dashed line, respectively), for the nonco- 
variant CCZ4 formulation with and without damping (black 
solid line and blue dotted line, respectively), and for the BSS- 
NOK formulation (red dashed line). We also report the dif- 
ferent values of coefficient / in the shift Eq. (fJTJ, which does 
change the growth rate of the unstable simulations, but does 
not remove the instability in the case of the fully covariant 
formulation 1 . The data refers to simulations having a coarse 
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Kj=l/2, k,=0.1/M, f=3/4 
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FIG. 4: L2-norm of the Hamiltonian constraint for the noncovariant 
CCZ4 formulation with and without damping terms (black solid line 
and blue dotted line, respectively), for the covariant CCZ4 formula- 
tion with and without damping terms (light-blue dot-dashed line and 
magenta long-dashed line, respectively), and for the BSSNOK for- 
mulation (red dashed line). Also indicated are the different values of 
coefficient / in the shift Eq. J2U . which however do not introduce 
qualitatively different behaviors. The data refers to the a simulations 
having a coarse resolution of ho/M = 0.48 and outer boundary at 
i? out = 2192.16 M. 



resolution of hg/M = 0.48 and outer boundary placed at 
-Rout = 2192.16 M, but similar behaviors have been seen also 
at higher and lower resolutions. 

Note that as the initial data settles and the evolution pro- 
ceeds, the CCZ4 formulation shows a violation of the Hamil- 
tonian constraint smaller than for the BSSNOK case (the L2- 
norm being at least 1 order of magnitude smaller), hence 
yielding a more accurate solution of the Einstein equations. 
However, after this initial stage, the evolutions with the CCZ4 
formulation can be considerably different according to the 
choice made for the parameters K3 and K\, More specifically, 
the covariant and damped system (i.e. K3 = 1, /ti ^ 0) ex- 
hibits a very rapid violation of the constraint at ~ 100 M and 
inevitably leads to a code crash (light-blue dot-dashed line 
in Fig. |4). Other variants of the CCZ4 formulation, on the 
other hand, show a different behavior. In particular, both of 
the undamped CCZ4 formulations (i.e. K3 = 1/2, 1, Ki = 0) 
lead to a successful merger, which can be easily identifiable 
as the peak at about ~ 350 — 380 M, and which is due to 
larger local violations of the constraints as the merger takes 



1 We have performed simulations also with K3 = \,K% = 0.1/M, / = 1, 
or ft 3 = 1, rei = 0.1/M, / = 3/4, and re 3 = 1, ki = 0, / = 3/4; in all 



cases we have found an instability (although with different growth rates), 
which we do not report to avoid overloading Fig. [4] 
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place 2 . At the same time, however, both implementations 
show a growth of the constraint violation (blue dotted line and 
magenta long-dashed line). This growth can be rather slow in 
the case / = 1, but it is likely to yield unstable evolutions 
on very long timescales. Finally, Fig. [4] shows that a nonco- 
variant and damped implementation of the CCZ4 formulation 
(i.e. «3 = 1/2, Ki 0; black solid line) leads not only to a 
stable merger and subsequent evolution, but it also provides a 
violation of the constraints which is at least 1 order of mag- 
nitude smaller than the corresponding one obtained with the 
BSSNOK evolution (red dashed line). This is one the main 
results of this paper and the ultimate justification for investi- 
gating this new formulation of the equations. 

We note that the behavior of the constraints described above 
for the CCZ4 formulation is indeed very similar to what al- 
ready experienced by many groups implementing the GH for- 
mulation 3 . In that case, in fact, the addition of the damp- 
ing terms was crucial to achieve stable black-hole evolu- 
tions id US EH. Altogether, the evolution shown in Fig. [Hal- 
ready provides the needed evidence that the new CCZ4 formu- 
lation, once suitable damping terms are added and the bound- 
ary conditions do not play a role, represents a considerable im- 
provement over the standard BSSNOK formulation. In what 
follows we will show that this continues to be the case also 
when the outer boundaries are chosen to produce incorrect 
data, or when they are placed very close to the merging bi- 
nary. 

Figure|5]reports with black solid lines the I = m = 2 mode 
of the gravitational waveform "J 4 extracted at r = 100 M for 
simulations having a coarse resolution ho/M = 0.60 and an 
outer boundary which is causally connected and at i? ou t = 
350.40 M (cf. Tabled. The top panel, in particular, refers to 
a simulation using the noncovariant and damped implemen- 
tation of the CCZ4 formulation (i.e. Z4d, with k 3 = 1/2, 
«i / 0), while the bottom one to a simulation using the 
BSSNOK formulation. Also shown with red dashed lines are 
the corresponding waveforms obtained when the outer bound- 
ary is causally disconnected and at i? ut = 2192.40 ill. As 
shown more clearly in the two insets, the CCZ4 formulation 
yields waveforms which are essentially identical and are un- 
affected by the constraint-violating outer boundaries. This is 
to be contrasted with the evolution performed with the BSS- 
NOK formulation and which shows strong signs of reflection 
att ~ 510 M. 

The reason behind this different behavior is to be found in 
the different way in which the two formulations handle the 
constraint- violations coming from the outer boundaries and is 
best appreciated in Fig. [6] where we show again the L2-norm 
of the ADM energy for the noncovariant and damped imple- 
mentation of the CCZ4 formulation (i.e. Z4d with K3 = 1/2, 
Ki ^ 0) and for the BSSNOK formulation. Note that both suf- 



2 Note that the time of merger is a gauge dependent quantity and can there- 
fore take place at slightly different times in different formulations. 

3 We recall that GH formulation can can be seen as a reduction of the Z4 
formalism 12711. 
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FIG. 5: Real part of the i = m = 2 mode of the gravitational wave- 
form ^4 (black solid line) extracted at r — 100 M for simulations 
having a coarse resolution ho/M = 0.60 and an outer boundary 
which is causally connected and at R ou t = 350.40 M. The top 
panel refers to a simulation using the noncovariant and damped im- 
plementation of the CCZ4 formulation (i.e. = 1/2, «i 0), 
while the bottom one to a simulation using the BSSNOK formula- 
tion; also shown are the corresponding waveforms obtained when 
flout = 2192.40 M (red dashed lines). 



fer of a very large increase at t ~ 250 M when the waves from 
the initial gauge settling of the binary, propagating at a speed 
of v g ~ v2, reach the outer boundary at i? out = 350.40 M 
and lead to larger violations. Also note that this increase in 
the constraint violation happens much earlier than the one as- 
sociates with the merger (which is at t ~ 350 M). As evi- 
dent from Fig.|6l the CCZ4 is able to recover efficiently from 
this violation, and the damping terms act in such a way that 
by t ~ 400 M the violation is completely removed, with the 
Hamiltonian constraint brought back to its minimum value. 
By contrast, the evolution with the BSSNOK formulation 
never recovers from the boundary contamination, leading to 
an increasing violation responsible for the incorrect behavior 
discussed in Fig. [5] The CZZ4 formulation experiences an- 
other increase in the violation at t ~ 750 M, when the gauge 
waves coming from the binary reach again the outer bound- 
ary, but once again the constraint damping terms act so as to 
remove the violation. 

An additional and concluding evidence of the constraint- 
damping properties of the CCZ4 formulation is shown is 
Fig. [7] where we report the evolution of the L2-norm of 
the Hamiltonian constraint (top panel) and of the root-mean- 
square of the momentum constraint (bottom panel) for the 
noncovariant and damped implementation of the CCZ4 for- 
mulation (i.e. Z4d with K3 = 1/2, K\ 7^ 0, black solid 
lines), and for the BSSNOK formulation (red dashed lines). 
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FIG. 6: L2-norm of the Hamiltonian constraint for the noncovariant 
and damped implementation of the CCZ4 formulation (i.e. Z4d with 
K3 = 1/2, Ki 0), and for the BSSNOK formulation (red dashed 
line). The data refers to the a simulations having a coarse resolution 
of ho/M = 0.60 and outer boundary placed at R ou t = 350.40 M. 



FIG. 7: L2-norm of the Hamiltonian constraint (top panel) and of the 
root-mean-square of the momentum constraint (bottom panel) for the 
noncovariant and damped implementation of the CCZ4 formulation 
(i.e. Z4d with k$ — 1/2, fti 7^ 0, black solid lines), and for the BSS- 
NOK formulation (red dashed lines). The data refers to simulations 
having a coarse resolution of ho/M — 1.20 and outer boundary at 
Rout = 199.20 M. 



The data refers to simulations performed with a plain Carte- 
sian outer boundary which is very close to the binary and at 
^out = 199.20 M (cf. Table|I]i. As in the previous figure, also 
here it is possible to detect the increase of the constraint vi- 
olations when gauge waves from the binary have reached the 
outer boundary at t ~ 140 M. 

Also in this case, the damping terms in the equations re- 
move rapidly the violations, which decay exponentially to 
their minimum values. Because the boundary is so close-in, 
this behavior of rapid increase and exponential decay takes 
place at least 3 times, both for the Hamiltonian and momen- 
tum constraints. Any formulation of the Einstein equations 
having this type of behavior is obviously preferable over one 
in which the violations are trapped in the computational do- 
main and are not allowed to be damped. 



IV. CONCLUSIONS 

By starting from the Z4 formulation l27ll and by includ- 
ing all the nonprincipal terms coming from the covariant form 
of the equations, we have introduced the CCZ4 formulation, 
i.e. the conformal and covariant formulation of the Z4 system, 
and proposed it as a new and effective way to solve numeri- 
cally the Einstein equations in arbitrary spacetimes. 

The new set of equations combines the most important fea- 
tures of the commonly used formulations of the Einstein equa- 
tions employed in numerical-relativity calculations. In partic- 
ular, it is able to make use of well-tested and robust gauge 
conditions which remove the need of excision and, at the same 



time, it is able to control dynamically the violation of the con- 
straint equations and to rapidly suppress them when they oc- 
cur. 

We have validated the robustness of the CCZ4 evolution 
system by performing a number of tests both in flat and in 
black-hole spacetimes. We have thus found that the CCZ4 
formulation without damping terms does not pass the stan- 
dard gauge-advection test, in analogy with the behavior of the 
BSSNOK formulation. However, when the damping terms are 
switched on, the new CCZ4 formulation passes the test stably 
and accurately. 

This ability of the formulation to control and damp vio- 
lations in the constraint equations has been confirmed also 
through the simulation of nonspinning black-hole binaries, 
which have been followed for about three orbits before merg- 
ing to a rapidly rotating black hole. Through a series of sim- 
ulations at different resolutions and with different treatments 
of the outer boundary - handled either with multiblocks and 
placed at a causally-disconnected distance, or with a Cartesian 
box and placed close to the binary - we have shown that not 
all of the implementations of the CCZ4 formulation lead to 
stable evolutions of binary black-hole spacetimes. 

Rather, we have found that the covariant form of the CCZ4 
formulation, in conjunction with the use of damping terms, 
leads to exponentially growing modes that rapidly destroy the 
numerical solution. Fortunately, the use of a noncovariant for- 
mulation and of damping terms leads not only to a stable evo- 
lution, but it also provides a violation of the constraints which 
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is at least 1 order of magnitude smaller than the corresponding 
one obtained with the BSSNOK evolution. A close compar- 
ison with simulations performed with the BSSNOK formula- 
tion using the same numerical setup, has also revealed that 
the CCZ4 formulation can efficiently recover from large vio- 
lations of the constraints, with the damping terms rapidly re- 
moving constraint violations produced at the outer boundary. 
By contrast, evolutions with the BSSNOK formulation expe- 
riencing similar violations never recover from the boundary 
contamination, leading to an increasing violation and incor- 
rect gravitational waves. 

Because the changes necessary to implement the new con- 
formal formulation in BSSNOK codes and the additional 
computational costs are very small, we propose the new for- 
mulation as a new standard for the numerical solution of the 
Einstein equations in generic 3D spacetimes. We expect, in 
fact, that a numerical solution of the Einstein equations hav- 
ing smaller violations of the constraints will also yield a more 
accurate modelling of the gravitational-wave emission, both 
in vacuum and nonvacuum spacetimes. 

At the same time, however, much remains to be done to 
fully understand the role played by the damping coefficients in 
fully nonlinear regimes and in the covariant form of the CCZ4 
formulation. Our experience with binary black-hole space- 
times has revealed, in fact, that there are situations in which 
the damping of the constraints interferes negatively with a 
fully covariant form of the CCZ4 formulation, leading to un- 



stable evolutions. In these cases, even small changes in the 
covariant character of the equations (e.g., by using k 3 = 0.9 
instead of K3 = 1) allows one to use nonzero damping co- 
efficients and hence to obtain a smaller violation of the con- 
straints. A systematic investigation of the space of parameters 
Ki x «2 x «3 is difficult due to the large computational costs 
of these simulations, but is clearly needed for a deeper under- 
standing of the behavior of the CCZ4 formulation. Much of 
our future work will be dedicated to elucidate this point. 
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